function [ Data ] = plotVIIVelocity( Data, Config, printPath )

%PLOTVIIVELOCITY Plot Vectrino II velocity data as histograms and mean profiles, returning mean 
% and median profiles in an updated Data structure.
%
%   Data = plotVIIVelocity( Data, Config, printPath) plots histograms and mean profiles of
%	velocity data for each velocity component with color coded range cells.
%
%	Inputs are the Data and Config structures generated during Matlab data export by the
%	Vectrino II acquistion software (also called MIDAS). An optional third argument
%	specifies a relative or absolute path for printing plots as PNGs.
%
%   See also plotVIISNR, plotVIIAmp, plotVIICorrelation.

%Copyright (c) <2012>, <Peter J. Rusello, pjrusello@gmail.com>
%All rights reserved.
%
%Redistribution and use in source and binary forms, with or without
%modification, are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright
%      notice, this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright
%      notice, this list of conditions and the following disclaimer in the
%      documentation and/or other materials provided with the distribution.
%    * Neither the name of the <organization> nor the
%      names of its contributors may be used to endorse or promote products
%      derived from this software without specific prior written permission.
%
%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
%ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
%WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
%DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER BE LIABLE FOR ANY
%DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
%(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
%LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
%ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
%SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

if nargin == 2
	printPath = '';
end
numberOfSamples = length( Data.Profiles_HostTime );

if isfield( Config, 'nCells' )
	nCells = Config.nCells;
	originalPrefix = '';
elseif isfield( Config, 'Original_nCells' )
	nCells = Config.Original_nCells;
	originalPrefix = 'Original_';
end

if Config.( [ originalPrefix 'coordSystem' ] ) == 1
	histogramBins = [ -1 * Config.( [ originalPrefix 'horizontalVelocityRange' ] ) / 1000:...
						0.001:...
						Config.( [ originalPrefix 'horizontalVelocityRange' ] ) / 1000 ];
else
	histogramBins = [ -1 * Config.( [ originalPrefix 'verticalVelocityRange' ] ) / 1000:...
						0.001:...
						Config.( [ originalPrefix 'verticalVelocityRange' ] ) / 1000 ];
end

histogramCounts = NaN * zeros( length( Data.Profiles_Range ), length( histogramBins ), 4 );
if Config.( [ originalPrefix 'coordSystem' ] ) == 1
	Data.Profiles_VelXmean = NaN * zeros( nCells, 3 );
	Data.Profiles_VelYmean = NaN * zeros( nCells, 3 );
	Data.Profiles_VelZ1mean = NaN * zeros( nCells, 3 );
	Data.Profiles_VelZ2mean = NaN * zeros( nCells, 3 );
else
	for beam = 1:4
		Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] ) = NaN * zeros( nCells, 3 );
	end
end

if isfield( Data, 'Profiles_DataQualityAllBeams' )
	if Config.( [ originalPrefix 'coordSystem' ] ) == 1
		Data.Profiles_VelXmean = NaN * zeros( nCells, 3, 2 );
		Data.Profiles_VelXmedian = NaN * zeros( nCells, 1 );
		Data.Profiles_VelYmean = NaN * zeros( nCells, 3, 2 );
		Data.Profiles_VelYmedian = NaN * zeros( nCells, 1 );
		Data.Profiles_VelZ1mean = NaN * zeros( nCells, 3, 2 );
		Data.Profiles_VelZ1median = NaN * zeros( nCells, 1 );
		Data.Profiles_VelZ2mean = NaN * zeros( nCells, 3, 2 );
		Data.Profiles_VelZ2median = NaN * zeros( nCells, 1 );
	else
		for beam = 1:4
			Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] ) = NaN * zeros( nCells, 3, 2 );
			Data.( [ 'Profiles_VelBeam' num2str( beam ) 'median' ] ) = NaN * zeros( nCells, 1 );
		end
	end
	histogramCountsScreened = NaN * zeros( length( Data.Profiles_Range ), length( histogramBins ), 4 );
end

velocityFieldNames = { 'X', 'Y', 'Z1', 'Z2' };
for beam = 1:4
    for cellNumber = 1:nCells
    	if Config.( [ originalPrefix 'coordSystem' ] ) == 1 % this is XYZ
	        Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 2, 1 ) = mean( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )( :, cellNumber ) );
	        Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'median' ] )( cellNumber, 2, 1 ) = median( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )( :, cellNumber ) );

	        if cellNumber == 1
		        [ bootstrapStatistics, bootstrapSamples ] = bootstrp( 1000, @mean, Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )( :, cellNumber ) );
		        bootstrapStatistics = sort( bootstrapStatistics );
		    else
				bootstrapStatistics = sort( mean( reshape( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )( bootstrapSamples, cellNumber ), numberOfSamples, 1000 ) ) );
			end
	        Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 1, 1 ) = bootstrapStatistics( 25 );
	        Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 3, 1 ) = bootstrapStatistics( 975 );

	        histogramCounts( cellNumber, :, beam ) = hist( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )( :, cellNumber ), histogramBins );
	        % if some type of screening has been performed affecting all velocity components, can also calculate statistics on the screened velocity data
	        if isfield( Data, 'Profiles_DataQualityAllBeams' )

		        histogramCountsScreened( cellNumber, :, beam ) = hist( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )...
		        	( Data.Profiles_DataQualityAllBeams( :, cellNumber ), cellNumber ), ...
		        	histogramBins );
		        Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 2, 2 ) = ...
		        	nanmean( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )...
		        	( Data.Profiles_DataQualityAllBeams( :, cellNumber ), cellNumber ) );
		        % feel free to uncomment the following lines to get bootstrap confidence
		        % intervals on the mean, be aware this is quite slow though since the
		        % bootsamples have to generated for each cell since screening is done on
		        % each individual cell
		        %bootstrapSamples = sort( bootstrp( 1000, @nanmean, Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )...
		        %	( Data.Profiles_DataQualityAllBeams( :, cellNumber ), cellNumber ) ) );
		        %Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 1, 2 ) = bootstrapSamples( 25 );
		        %Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 3, 2 ) = bootstrapSamples( 975 );

		    end
    	else % this is Beam
	        Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 2, 1 ) = mean( Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )( :, cellNumber ) );
	        Data.( [ 'Profiles_VelBeam' num2str( beam ) 'median' ] )( cellNumber, 1 ) = median( Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )( :, cellNumber ) );
	        if cellNumber == 1
		        [ bootstrapStatistics, bootstrapSamples ] = bootstrp( 1000, @mean, Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )( :, cellNumber ) );
		        bootstrapStatistics = sort( bootstrapStatistics );
		    else
				bootstrapStatistics = sort( mean( reshape( Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )( bootstrapSamples, cellNumber ), numberOfSamples, 1000 ) ) );
			end
	        Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 1, 1 ) = bootstrapStatistics( 25 );
	        Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 3, 1 ) = bootstrapStatistics( 975 );

	        histogramCounts( cellNumber, :, beam ) = hist( Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )( :, cellNumber ), histogramBins );
	        % if some type of screening has been performed affecting all velocity components, can also calculate statistics on the screened velocity data
	        if isfield( Data, 'Profiles_DataQualityAllBeams' )

		        histogramCountsScreened( cellNumber, :, beam ) = hist( ...
		        	Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )...
		        		( Data.( [ 'Profiles_DataQualityBeam' num2str( beam ) ] )( :, cellNumber ), cellNumber ), ...
		        		histogramBins );
		
		        Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 2, 2 ) = nanmean( ...
		        	Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )...
		        	( Data.( [ 'Profiles_DataQualityBeam' num2str( beam ) ] )( :, cellNumber ), cellNumber ) );
		        Data.( [ 'Profiles_VelBeam' num2str( beam ) 'median' ] )( cellNumber, 2 ) = nanmedian( ...
		        	Data.( [ 'Profiles_VelBeam' num2str( beam ) ] )...
		        	( Data.( [ 'Profiles_DataQualityBeam' num2str( beam ) ] )( :, cellNumber ), cellNumber ) );
		        % feel free to uncomment the following lines to get bootstrap confidence
		        % intervals on the mean, be aware this is quite slow though since the
		        % bootsamples have to generated for each cell since screening is done on
		        % each individual cell
		        %bootstrapSamples = sort( bootstrp( 1000, @nanmean, ...
		        %	Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } ] )...
		        %	( Data.( [ 'Profiles_DataQualityBeam' num2str( beam ) ] )( :, cellNumber ), cellNumber ) ) );
		        %Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 1, 2 ) = bootstrapSamples( 25 );
		        %Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 3, 2 ) = bootstrapSamples( 975 );

		    end

	    end        
    end
end


figure
cellColors = colormap( jet( nCells ) );

for beam = 1:4
    subplot( 1, 4, beam )
    for cellNumber = 1:nCells
    	if Config.( [ originalPrefix 'coordSystem' ] ) == 1
	        plot( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 2, 1 ), Data.Profiles_Range( cellNumber ), '*', 'color', cellColors( cellNumber, : ) )
    	    hold on
    	else
    		plot( Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 2, 1 ), Data.Profiles_Range( cellNumber ), '*', 'color', cellColors( cellNumber, : ) )
    	    hold on
		end

		if isfield( Data, 'Profiles_DataQualityAllBeams' )
			if Config.( [ originalPrefix 'coordSystem' ] ) == 1
				plot( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( cellNumber, 2, 2 ), Data.Profiles_Range( cellNumber ), 'o', 'color', cellColors( cellNumber, : ) )
				hold on
			else
				plot( Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( cellNumber, 2, 2 ), Data.Profiles_Range( cellNumber ), 'o', 'color', cellColors( cellNumber, : ) )
				hold on
			end
		end
		
    end
    
    if Config.( [ originalPrefix 'coordSystem' ] ) == 1
		plot( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( :, 1, 1 ), Data.Profiles_Range, 'k--' )
		plot( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( :, 3, 1 ), Data.Profiles_Range, 'k--' )
	else
		plot( Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( :, 1, 1 ), Data.Profiles_Range, 'k--' )
		plot( Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( :, 3, 1 ), Data.Profiles_Range, 'k--' )
	end	
	
    if isfield( Data, 'Profiles_DataQualityAllBeams' )
		if Config.( [ originalPrefix 'coordSystem' ] ) == 1
			plot( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( :, 1, 2 ), Data.Profiles_Range, 'k--' )
			plot( Data.( [ 'Profiles_Vel' velocityFieldNames{ beam } 'mean' ] )( :, 3, 2 ), Data.Profiles_Range, 'k--' )
		else
			plot( Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( :, 1, 2 ), Data.Profiles_Range, 'k--' )
			plot( Data.( [ 'Profiles_VelBeam' num2str( beam ) 'mean' ] )( :, 3, 2 ), Data.Profiles_Range, 'k--' )
		end	
	end
    grid on
    set( gca, 'ydir', 'reverse' )
    if Config.( [ originalPrefix 'coordSystem' ] ) == 1
	    xlabel( [ 'Vel. ' velocityFieldNames{ beam } ' (m/s)' ] )
	else
	    xlabel( [ 'Vel. B_' num2str( beam ) ' (m/s)' ] )
	end
    if beam == 1
        ylabel( 'Range (m)' )
    end
    set( gca, 'ylim', [ Data.Profiles_Range( 1 ) Data.Profiles_Range( end ) ] )
end

if ~isempty( printPath )
	if Config.( [ originalPrefix 'coordSystem' ] ) == 1
		print( '-dpng', '-r300', [ printPath 'meanXYZVelocityProfile.png' ] )
	else
		print( '-dpng', '-r300', [ printPath 'meanBeamVelocityProfile.png' ] )
	end
end

figure
for beam = 1:4
    subplot( 4, 1, beam )
    for cellNumber = 1:nCells
        nonzeroBins = histogramCounts( cellNumber, :, beam ) > 0;
        plot( histogramBins( nonzeroBins ), histogramCounts( cellNumber, nonzeroBins, beam ) / numberOfSamples, 'color', cellColors( cellNumber, : ) )
        hold on
        if cellNumber == 1
            minBin( beam ) = min( histogramBins( nonzeroBins ) );
            maxBin( beam ) = max( histogramBins( nonzeroBins ) );
        end
        if min( histogramBins( nonzeroBins ) ) < minBin( beam )
            minBin( beam ) = min( histogramBins( nonzeroBins ) );
        end
        if max( histogramBins( nonzeroBins ) ) > maxBin( beam )
            maxBin( beam ) = max( histogramBins( nonzeroBins ) );
        end
    end
    grid on
    set( gca, 'xlim', [ minBin( beam ) maxBin( beam ) ] )
    if Config.( [ originalPrefix 'coordSystem' ] ) == 1 & beam == 4
	    xlabel( [ 'Vel. ' velocityFieldNames{ beam } ' (m/s)' ] )
	elseif beam == 4
	    xlabel( [ 'Vel. B_' num2str( beam ) ' (m/s)' ] )
	end
    if beam == 3
        ylabel( 'Occurrence (%)' )
    end

end


if ~isempty( printPath )
	if Config.( [ originalPrefix 'coordSystem' ] ) == 1
		print( '-dpng', '-r300', [ printPath 'XYZVelocityHistogram.png' ] )
	else
		print( '-dpng', '-r300', [ printPath 'beamVelocityHistogram.png' ] )
	end
end
